Análisis de Series Temporales

Introducción

Las series temporales son técnicas de aprendizaje que permiten analizar los datos y la secuencia de tiempo para predecir eventos futuros. Permiten explicar la evolución histórica de una variable a lo largo del tiempo y predecir sus valores futuros. En este trabajo se presentan un conjunto de modelos basados en análisis temporales que permiten modelar y predecir el comportamiento de la contratación de hipotecas en España, específicamente en el País Vasco. Para ello se cuenta con un dataset que contiene datos de la cantidad de hipotecas contratadas en cada provincia del País Vasco desde 2003 hasta 2019. El objetivo es desarrollar un modelo de series temporales sobre los datos de 2003 a 2018 que minimice el error de predicción en los datos de 2019.

Desarrollo

En esta sección se construye la serie correspondiente a la comunidad autónoma. Asimismo, se presentan y analizan los diferentes modelos ARMA de series temporales, mostrando los modelos intermedios y los razonamientos que condujeron a proponer cada uno de ellos.

Análisis exploratorio.

Inicialmente se almacenan en memoria los datos correspondientes a la comunidad autónoma País Vasco. El conjunto de datos original contiene la cantidad de hipotecas contratadas mensualmente en cada provincia, con el objetivo de realizar el análisis a nivel de comunidad, se suman las cantidades de cada provincia obteniendo una única columna total que representa la cantidad de contrataciones mensuales en la comunidad.

original <- readxl::read_xlsx("./data/Hipotecas.xlsx")
original <- original[c('...1', '01 Álava', '20 Gipuzkoa', '48 Vizcaya')]
colnames(original) <- c('date', 'Alava', 'Gipuzkoa', 'Vizcaya')
original$date <- as.Date(paste(original$date, "01", sep = "-"), "%YM%m-%d")
original$total <- original$Alava + original$Gipuzkoa + original$Vizcaya
data <- original[c('date', 'total')]

Para realizar una validación correcta de los modelos a desarrollar el conjunto de datos se divide en entrenamiento y test. Los datos desde 2003 hasta 2018 se utilizarán para entrenar los modelos, mientras que los datos de 2019 serán utilizados en la etapa de validación.

data.train <- subset(data, year(data$date) != 2019)
data.test <- subset(data, year(data$date) == 2019)
data.train.ts <- as.ts(data.train$total, frequency = 12)
data.test.ts <- as.ts(data.test$total, frequency=12)

En la siguiente gráfica se muestra la serie temporal que analizaremos en este trabajo, obtenida tras la suma de los datos de las tres provincias que componen la comunidad.

ggplot(aes(x= date, y = total), data = data.train) +
  geom_line(color = '#d84519', size = 1) + 
  xlab('FECHA') + ylab('Hipotecas')

A simple vista se puede determinar que en la serie hay presencia de datos atípicos tipo pulso, por ejemplo, el pico observable en 2006, que rompe de manera puntual el patrón de la serie. Asimismo, se observa que la serie no es estacionaria pues se evidencia que no es, al menos, estacionaria en media. Comprobemos las hipótesis de estacionalidad de manera analítica en la siguiente sección.

Modelos

El objetivo es lograr que la serie sea estacionaria para que todos los instantes sean comparables. Inicialmente analizaremos si la serie es estacionaria en varianza. Se evalúa la necesidad de transformar la serie para hacerla estacionaria en varianza a través de la transformación BoxCox.

box_cox <- boxcox(total ~ date,
                  data = data.train,
                  lambda = c(0, 0.5, 1))

lambda <- box_cox$x[which.max(box_cox$y)]
lambda

## [1] 0.6464646

Al obtener un lambda de 0.65, diferente de 1, indica que la serie no es estacionaria en varianza, por lo tanto, hay que aplicar la transformación elevando la serie a lambda.

data.train$exp_total = data.train$total**lambda
data.test$exp_total = data.test$total**lambda
data.train.ts <- as.ts(data.train$exp_total, frequency = 12)
data.test.ts <- as.ts(data.test$exp_total, frequency=12)

Tras esta transformación ya podemos afirmar que la serie es estacionaria en varianza. Verifiquemos si el modelo es estacionario en media a través del test Dickey-Fuller.

adf.test(data.train.ts)

## 
##  Augmented Dickey-Fuller Test
## 
## data:  data.train.ts
## Dickey-Fuller = -1.9115, Lag order = 5, p-value = 0.6132
## alternative hypothesis: stationary

Al aplicar el test de Dickey-Fuller se obtiene un p-value alto, por lo que no se puede rechazar la hipótesis nula que plantea que la serie no es estacionaria en media. Esto sugiere que podría ser necesaria una diferencia regular. En este caso, ya que el test aplicado tiene una baja potencia se decidirá más adelante la necesidad de aplicar o no una diferencia.

Ajuste de modelos

El objetivo que se persigue es identificar el proceso que subyace bajo los datos, lo cual consiste en identificar los órdenes p y q del modelo ARMA que generó la serie temporal. Esto se puede determinar a través de las funciones de autocorrelación simple (f.a.s) y parcial (f.a.p).

A continuación, se visualizan la funciones de autocorrelación simple y parcial de residuos para comprobar si se cumple la hipótesis de ruido blanco, que implica que los residuos están incorrelados entre sí. En este caso, se construye sobre los datos originales porque la serie original inicialmente es toda residuo.

En la siguiente gráfica se muestra la función de correlación simple f.a.s. Se puede notar que la cantidad de palitos que están fuera de la banda son infinitos por lo que estamos en presencia de un MA (denota el orden de la q) infinito. En clase fue demostrado que un MA infinito es equivalente a un AR (parte autorregresiva, denota el orden de la p) finito, por lo tanto es necesario determinar el orden del AR finito.

acf(data.train.ts, lag.max = 48, xlab = "Hipotecas",
    main= "Función de autocorrelacion simple")

Para determinar el orden de la AR se dibuja el gráfico de autocorrelación parcial. De igual manera, el orden de la p se determina teniendo en cuenta la cantidad de líneas fuera de la banda en la función f.a.p.

pacf(data.train.ts, lag.max = 48, xlab = "Hipotecas",
     main = "Funcion de autocorrelacion parcial")

La imagen sugiere, de momento, que el orden de la parte autorregresiva puede ser 1 porque es la línea que más sobresale sobre el resto. Por esta razón, se propone como modelo inicial un AR(1).

El gráfico de la función de autocorrelación simple f.a.s, aunque con un orden infinito, muestra también que existe cierta periodicidad en los múltiplos de 4 lo que sugiere que el período pudiera ser 4. Sin embargo, al tratarse de una serie mensual lo más natural es establecer período 12, un período 12 siempre reflejará también el período 4. Además, por cuestiones de explicabilidad de los modelos, seria más conveniente seleccionar período 12 puesto que es más evidente la relación que existe entre un mes determinado y el mismo mes del año anterior, que entre un mes y 4 meses anteriores.

Modelo 1: ARIMA(1, 0, 0)

fit.1 <- Arima(data.train.ts,
               order = c(1, 0, 0),
               seasonal = list(order=c(0,0,0), period = 12),
               method = 'ML')

fit.1

## Series: data.train.ts 
## ARIMA(1,0,0) with non-zero mean 
## 
## Coefficients:
##          ar1      mean
##       0.9013  139.9555
## s.e.  0.0306   12.4928
## 
## sigma^2 estimated as 322.8:  log likelihood=-826.87
## AIC=1659.75   AICc=1659.88   BIC=1669.52

coeftest(fit.1)

## 
## z test of coefficients:
## 
##             Estimate Std. Error z value  Pr(>|z|)    
## ar1         0.901331   0.030554  29.500 < 2.2e-16 ***
## intercept 139.955522  12.492826  11.203 < 2.2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Al analizar los coeficientes del modelo, se puede notar que los dos parámetros son significativos y distintos de 0. La suma de los parámetros no es cercana a 1 por lo que no se viola la condición de estacionalidad. Analicemos las correlaciones entre los parámetros para evitar que exista colinealidad.

cor.arma(fit.1)

##                   ar1   intercept
## ar1       1.000000000 0.006460326
## intercept 0.006460326 1.000000000

La correlación entre ambos parámetros no es alta por lo cual no existe colinealidad. A continuación, se analiza la hipótesis de ruido blanco. El siguiente test muestra si las correlaciones de los 6 primeros residuos son ruido blanco. Los valores muestran que aún no hemos logrado ruido blanco.

Box.test.2(residuals(fit.1),
           nlag = c(6, 12, 18, 24, 30, 36, 48))

##      Retard p-value
## [1,]      6       0
## [2,]     12       0
## [3,]     18       0
## [4,]     24       0
## [5,]     30       0
## [6,]     36       0
## [7,]     48       0

Visualicemos la afirmación anterior a través de las gráficas de autocorrelación. Efectivamente, en la siguiente figura, se puede observar que la serie no cumple la hipótesis de ruido blanco pues existen múltiples líneas fuera de la banda.

acf(fit.1$residuals, lag.max = 100, xlab = "Hipotecas",
    main= "Función de autocorrelacion simple")

Se puede observar cierta simetría lo que indica que hay una parte regular que va asociada a los instantes inmediatamente anteriores y una parte estacional que va asociada a los mismos instantes pero retardados ese período de tiempo. Esto sugiere la existencia de un modelo SARIMA multiplicativo.

Continuamos teniendo un MA infinito, por lo tanto seguimos teniendo un AR finito. Para determinar el orden de la p dibujamos nuevamente en gráfico de autocorrelación parcial.

pacf(fit.1$residuals, lag.max = 100, xlab = "Hipotecas",
     main = "Funcion de autocorrelacion parcial")

En el gráfico de la función f.a.p el palito que más sobresale por fuera de la banda es el 1 lo cual sugierela inclusión de un modelo AR(1) en la parte estacional, que habíamos decidido adicionar con anterioridad.

Modelo 2: ARIMA(1, 0, 0)x(1, 0, 0)[12]

fit.2 <- Arima(data.train.ts,
               order = c(1, 0, 0),
               seasonal = list(order = c(1, 0, 0), period=12),
               method = 'ML')

fit.2

## Series: data.train.ts 
## ARIMA(1,0,0)(1,0,0)[12] with non-zero mean 
## 
## Coefficients:
##          ar1    sar1      mean
##       0.7739  0.5778  141.1561
## s.e.  0.0548  0.0689   10.3325
## 
## sigma^2 estimated as 228.8:  log likelihood=-795.4
## AIC=1598.81   AICc=1599.02   BIC=1611.84

coeftest(fit.2)

## 
## z test of coefficients:
## 
##             Estimate Std. Error z value  Pr(>|z|)    
## ar1         0.773899   0.054816 14.1182 < 2.2e-16 ***
## sar1        0.577813   0.068931  8.3825 < 2.2e-16 ***
## intercept 141.156064  10.332487 13.6614 < 2.2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Los parámetros son significativos y en ninguno de los casos el 1 pertenece al intervalo de confianza, por lo tanto, se cumple la condición de estacionalidad. Analicemos la correlación entre parámetros.

cor.arma(fit.2)

##                   ar1        sar1   intercept
## ar1        1.00000000 -0.57640143 -0.00175156
## sar1      -0.57640143  1.00000000  0.02249438
## intercept -0.00175156  0.02249438  1.00000000

Al analizar las correlaciones, se puede notar que no existe una correlación significativamente alta entre los parámetros que puede afectar el modelo. Analicemos la existencia de ruido blanco en el modelo a través de las correlaciones y la función f.a.s.

Box.test.2(residuals(fit.2),
           nlag = c(6, 12, 18, 24, 30, 36, 48))

##      Retard    p-value
## [1,]      6 0.00002063
## [2,]     12 0.00000387
## [3,]     18 0.00001581
## [4,]     24 0.00001803
## [5,]     30 0.00008436
## [6,]     36 0.00012132
## [7,]     48 0.00071927

acf(fit.2$residuals, lag.max = 48, xlab = "Hipotecas",
     main = "Funcion de autocorrelacion simple")

Las autocorrelaciones indican que el p-value de los 6 primeros son muy cercanos a cero, por lo tanto, no existe aún ruido blanco. Esto se reafirma si analizamos la gráfica de autocorrelación simple. La máxima autocorrelación simple distinta de 0 determina el orden de MA, por eso, en este caso, se puede sugerir un modelo MA(1) en la parte regular. Como se trata de un MA finito no es necesario analizar la autocorrelación parcial y pasamos directamente a un nuevo modelo.

Modelo 3: ARIMA(1,0,1)(1,0,0)[12]

fit.3 <- Arima(data.train.ts,
               order = c(1, 0, 1),
               seasonal = list(order = c(1, 0, 0), period=12),
               method = 'ML')

fit.3

## Series: data.train.ts 
## ARIMA(1,0,1)(1,0,0)[12] with non-zero mean 
## 
## Coefficients:
##          ar1      ma1    sar1      mean
##       0.9803  -0.6603  0.4779  142.8461
## s.e.  0.0135   0.0544  0.0633   24.1849
## 
## sigma^2 estimated as 169.2:  log likelihood=-765.89
## AIC=1541.78   AICc=1542.11   BIC=1558.07

coeftest(fit.3)

## 
## z test of coefficients:
## 
##             Estimate Std. Error  z value  Pr(>|z|)    
## ar1         0.980294   0.013478  72.7328 < 2.2e-16 ***
## ma1        -0.660289   0.054351 -12.1485 < 2.2e-16 ***
## sar1        0.477922   0.063268   7.5540 4.222e-14 ***
## intercept 142.846076  24.184913   5.9064 3.496e-09 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

En este caso, todos los parámetros son significativos, sin embargo, el 1 pertenece al intervalo de confianza del parámetro asociado a la AR(1) lo que indica que no se está cumpliendo la condición de estacionalidad y es necesario transformar el AR(1) en una diferencia. De todas formas, analicemos la correlación y la hipótesis de ruido blanco.

cor.arma(fit.3)

##                   ar1         ma1         sar1    intercept
## ar1        1.00000000 -0.42927390 -0.273345876  0.047392350
## ma1       -0.42927390  1.00000000  0.109105050 -0.017882559
## sar1      -0.27334588  0.10910505  1.000000000 -0.001957786
## intercept  0.04739235 -0.01788256 -0.001957786  1.000000000

No existe una alta correlación entre los parámetros.

Box.test.2(residuals(fit.3),
           nlag = c(6, 12, 18, 24, 30, 36, 48))

##      Retard   p-value
## [1,]      6 0.7181104
## [2,]     12 0.5739264
## [3,]     18 0.6931901
## [4,]     24 0.5070736
## [5,]     30 0.6844058
## [6,]     36 0.4277117
## [7,]     48 0.5098180

acf(fit.3$residuals, lag.max = 48, xlab = "Hipotecas",
     main = "Funcion de autocorrelacion simple")

pacf(fit.3$residuals, lag.max = 48, xlab = "Hipotecas",
     main = "Funcion de autocorrelacion parcial")

Al analizar las autocorrelaciones se puede notar que la inclusión de un modelo MA(1) nos condujo a la presencia de ruido blanco, como muestra también la gráfica de autocorrelación simple. Apliquemos una diferencia en la parte regular, como habíamos determinado anteriormente, para que se cumpla la hipótesis de estacionalidad.

Modelo 4: ARIMA(0,1,1)(1,0,0)[12]

fit.4 <- Arima(data.train.ts,
               order = c(0, 1, 1),
               seasonal = list(order = c(1, 0, 0), period=12),
               method = 'ML')
fit.4

## Series: data.train.ts 
## ARIMA(0,1,1)(1,0,0)[12] 
## 
## Coefficients:
##           ma1    sar1
##       -0.6742  0.4667
## s.e.   0.0505  0.0620
## 
## sigma^2 estimated as 169.2:  log likelihood=-761.82
## AIC=1529.63   AICc=1529.76   BIC=1539.39

coeftest(fit.4)

## 
## z test of coefficients:
## 
##       Estimate Std. Error  z value  Pr(>|z|)    
## ma1  -0.674160   0.050459 -13.3604 < 2.2e-16 ***
## sar1  0.466659   0.062031   7.5229 5.356e-14 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Al aplicar la diferencia, se ha eliminado la constante. Aunque en la mayoría de los casos al aplicar una diferencia este parámetro se hace no significativo, comprobemos si realmente la constante no es significtiva forzando su inclusión en el modelo.

fit.4 <- Arima(data.train.ts,
               include.constant = T,
               order = c(0, 1, 1),
               seasonal = list(order = c(1, 0, 0), period=12),
               method = 'ML')
coeftest(fit.4)

## 
## z test of coefficients:
## 
##        Estimate Std. Error  z value  Pr(>|z|)    
## ma1   -0.676555   0.050596 -13.3717 < 2.2e-16 ***
## sar1   0.463351   0.062353   7.4311 1.077e-13 ***
## drift -0.300919   0.541296  -0.5559    0.5783    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Al inlcuir la constante, el p-valor de la misma es 0.5 por lo que efectivamente se puede determinar que no es significativa y lo podemos eliminar nuevamente.

fit.4 <- Arima(data.train.ts,
               order = c(0, 1, 1),
               seasonal = list(order = c(1, 0, 0), period=12),
               method = 'ML')
fit.4

## Series: data.train.ts 
## ARIMA(0,1,1)(1,0,0)[12] 
## 
## Coefficients:
##           ma1    sar1
##       -0.6742  0.4667
## s.e.   0.0505  0.0620
## 
## sigma^2 estimated as 169.2:  log likelihood=-761.82
## AIC=1529.63   AICc=1529.76   BIC=1539.39

coeftest(fit.4)

## 
## z test of coefficients:
## 
##       Estimate Std. Error  z value  Pr(>|z|)    
## ma1  -0.674160   0.050459 -13.3604 < 2.2e-16 ***
## sar1  0.466659   0.062031   7.5229 5.356e-14 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Todos los parámetros son significativos y la suma no es cercana a 1, por lo tanto, se vuelve a cumplir la hipótesis de estacionalidad. En el caso de MA(1) también se cumple la hipótesis de invertibilidad que aunque no es necesaria para ajustar un modelo ARMA sí es deseable pues su incumplimiento lleva a que datos más alejados tengan más peso en la predicción que datos más recientes. Analicemos las correlaciones en los parámetros para evitar la colinealidad en el modelo.

cor.arma(fit.4)

##             ma1       sar1
## ma1  1.00000000 0.06378629
## sar1 0.06378629 1.00000000

La matriz de correlación muestra que no hay peligro de autocolinealidad dada la baja correlación entre los parámetros. Analicemos si este nuevo modelo cumple la hipótesis de ruido blanco.

Box.test.2(residuals(fit.4),
           nlag = c(6, 12, 18, 24, 30, 36, 48))

##      Retard   p-value
## [1,]      6 0.7270230
## [2,]     12 0.6153085
## [3,]     18 0.6866355
## [4,]     24 0.4826159
## [5,]     30 0.6498918
## [6,]     36 0.3884261
## [7,]     48 0.5178972

acf(fit.4$residuals, lag.max = 48, xlab = "Hipotecas",
     main = "Funcion de autocorrelacion simple")

Al analizar el gráfico de autocorrelación simple y analizar las correlaciones se puede decir que ya estamos en presencia de ruido blanco. Sin embargo, aún existen unos pocos palitos que se salen de las bandas como el 19, 32, entre otros. Como no se observa un criterio evidente que permita corregir este hecho pasaremos a incluir variables explicativas.

Análisis de intervenciones

En una serie mensual hay tres efectos importantes que reflejar. El efecto calendario que consiste en reflejar en el modelo la diferencia entre los días festivos y laborables que hay en un mes. Además, se deben tener en cuenta los días de Semana Santa y los años bisiestos.

La siguiente función, tomada del ejemplo visto en clase, permite determinar la cantidad de festivos en un mes, cuántos días de Semana Santa hubieron en ese mes y si el mes corresponde a un año bisiesto o no.

calculoExplicativasCalendario <- function(variableFecha, domingoYFestivosJuntos){
  if (month(max(variableFecha)) %in% c(1,3,5,7,8,10,12)) {
    diasHastaFinMes <- 30
  } else if (month(max(variableFecha)) %in% c(4,6,9,11)) {
    diasHastaFinMes <- 29
  } else if (year(max(variableFecha))%%4==0) {
    diasHastaFinMes <- 28
  } else {diasHastaFinMes <- 27}
  
  todasLasFechas <- data.frame(fechas=seq(min(variableFecha),
                                          max(variableFecha)+diasHastaFinMes,
                                          by="days"))
  
  domingoResurrecion <- as.Date(Easter(year(min(variableFecha)):year(max(variableFecha))))
  lunesPascua <- domingoResurrecion+1
  sabadoSanto <- domingoResurrecion-1
  viernesSanto <- domingoResurrecion-2
  juevesSanto <- domingoResurrecion-3
  
  semanaSanta <- sort(c(juevesSanto, viernesSanto, sabadoSanto, domingoResurrecion, lunesPascua))
  
  semanaSanta <- data.frame(fechas=semanaSanta, semanaSanta=rep(1,length(semanaSanta)))
  
  todasLasFechas_2 <- merge(x = todasLasFechas, y = semanaSanta, by = "fechas", all.x = TRUE)
  
  todasLasFechas_2$semanaSanta[is.na(todasLasFechas_2$semanaSanta)] <- 0
  
  calendario <- todasLasFechas
  
  calendario$diaSemana <- as.factor(wday(calendario$fecha))
  calendario$diaMes <- as.factor(day(calendario$fecha))
  calendario$mes <- as.factor(month(calendario$fecha))
  calendario$anyo <- as.factor(year(calendario$fecha))
  
  calendario$p_01ene <- ifelse(calendario$diaMes==1 & calendario$mes==1, 1, 0)
  calendario$p_06ene <- ifelse(calendario$diaMes==6 & calendario$mes==1, 1, 0)
  calendario$p_19mar <- ifelse(calendario$diaMes==19 & calendario$mes==3, 1, 0)
  calendario$p_01may <- ifelse(calendario$diaMes==1 & calendario$mes==5, 1, 0)
  calendario$p_15ago <- ifelse(calendario$diaMes==15 & calendario$mes==8, 1, 0)
  calendario$p_12oct <- ifelse(calendario$diaMes==12 & calendario$mes==10,1, 0)
  calendario$p_01nov <- ifelse(calendario$diaMes==1 & calendario$mes==11, 1 ,0)
  calendario$p_06dic <- ifelse(calendario$diaMes==6 & calendario$mes==12, 1 ,0)
  calendario$p_08dic <- ifelse(calendario$diaMes==8 & calendario$mes==12, 1 ,0)
  calendario$p_25dic <- ifelse(calendario$diaMes==25 & calendario$mes==12, 1 ,0)
  
  calendario$festivo <- rowSums(subset(calendario, select=p_01ene:p_25dic))
  
  if (domingoYFestivosJuntos==0){
    
    calendario$sabado <- ifelse(calendario$diaSemana==7, 1 ,0)
    calendario$domingo <- ifelse(calendario$diaSemana==1, 1 ,0)
    
    calendario$laborable <- 1-calendario$sabado-calendario$domingo
    
  } else {
    
    calendario$sabado <- ifelse(calendario$diaSemana==7, 1 ,0)
    calendario$domingo <- ifelse(calendario$diaSemana==1, 1 ,0)
    calendario$domingo <- ifelse(calendario$domingo==1 | calendario$festivo==1, 1 ,0)
    
    calendario$laborable <- 1-calendario$sabado-calendario$domingo    
  }
  
  calendario_2 <- calendario[, c("fechas", "mes", "anyo", "sabado", "domingo", "laborable", "festivo")]
  
  todasLasFechasFinal <- merge(x = todasLasFechas_2, y = calendario_2,
                               by = "fechas", all.x = TRUE)
  
  calendarioAnyoMes <- aggregate(todasLasFechasFinal[,c("sabado","domingo",
                                                        "laborable", "semanaSanta", "festivo")],
                                 by=list(mes=todasLasFechasFinal$mes,
                                         anyo=todasLasFechasFinal$anyo),
                                 "sum")
  
  calendarioAnyoMes$dt <- calendarioAnyoMes$laborable-(5/2)*(calendarioAnyoMes$sabado+calendarioAnyoMes$domingo)
  
  calendarioAnyoMes$anyoNum <- as.numeric(levels(calendarioAnyoMes$anyo))[calendarioAnyoMes$anyo]
  
  calendarioAnyoMes$bisiesto <- ifelse(calendarioAnyoMes$mes==2 &(calendarioAnyoMes$anyoNum %% 4)==0, 1 ,0)
  
  if (domingoYFestivosJuntos==0){
    explicativasCalendario <- cbind(fecha=variableFecha, calendarioAnyoMes[, c("semanaSanta", "dt", "bisiesto", "festivo")])
  } else {
    explicativasCalendario <- cbind(fecha=variableFecha, calendarioAnyoMes[, c("semanaSanta", "dt", "bisiesto")])
  }
  
  return(explicativasCalendario)
  
}

La variable dt representa la variable calendario, si toma valores positivos indica el número de días festivos de más que tiene un mes respecto a un mes perfecto (como febrero) y si es negativa, el número días de menos respecto a un mes perfecto. Intentemos solamente con las variables calendario, los días de Semanas Santa y si el año es bisiesto o no.

exp.calendar.train <- calculoExplicativasCalendario(data.train$date,domingoYFestivosJuntos=0)

calendar.train <- 
  as.matrix(
    exp.calendar.train[,c("semanaSanta", "dt", "bisiesto")])

Se toma el último modelo que habíamos llegado a la conclusión que lograba ruido blanco y le incluimos los factores explicativos referentes al calendaria, Semana Santa y años bisiestos. A este modelo le llamaremos ARIMA(0,1,1)(1,0,0)[12] + calendario, aunque inclye los otros factores mencionados.

Modelo 5: ARIMA(0,1,1)(1,0,0)[12] + calendario

fit.calendar <- Arima(data.train.ts,
               order = c(0, 1, 1),
               seasonal = list(order = c(1, 0, 0), period=12),
               xreg = calendar.train,
               method = 'ML')
fit.calendar

## Series: data.train.ts 
## Regression with ARIMA(0,1,1)(1,0,0)[12] errors 
## 
## Coefficients:
##           ma1    sar1  semanaSanta      dt  bisiesto
##       -0.5970  0.4755      -2.3426  1.3606    2.6405
## s.e.   0.0599  0.0621       0.5776  0.2716    4.7327
## 
## sigma^2 estimated as 140.3:  log likelihood=-742.34
## AIC=1496.68   AICc=1497.14   BIC=1516.2

coeftest(fit.calendar)

## 
## z test of coefficients:
## 
##              Estimate Std. Error z value  Pr(>|z|)    
## ma1         -0.597044   0.059922 -9.9637 < 2.2e-16 ***
## sar1         0.475524   0.062067  7.6614 1.839e-14 ***
## semanaSanta -2.342564   0.577635 -4.0554 5.004e-05 ***
## dt           1.360636   0.271571  5.0102 5.436e-07 ***
## bisiesto     2.640477   4.732710  0.5579    0.5769    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

El hecho de que la variable calendario (dt) sea positiva indica, como es de esperar, que cuanto más días laborables más contratos de hipotecas. Mientras más días de Semana Santa haya en un mes, se realizan menos contrataciones. Todas las variables son significativas excepto la relativa a los años bisiestos, sin embargo, al ser positiva no afecta mucho al modelo y tiene sentido incluirla, porque si un año tiene un días más, es un día más donde se pueden comprar hipotecas.

Analicemos la correlación entre los parámetros.

cor.arma(fit.calendar)

##                     ma1        sar1 semanaSanta           dt     bisiesto
## ma1          1.00000000  0.01521880 -0.05592220  0.019007244  0.014101586
## sar1         0.01521880  1.00000000  0.01547795 -0.060104239 -0.056010132
## semanaSanta -0.05592220  0.01547795  1.00000000  0.074654784  0.029814131
## dt           0.01900724 -0.06010424  0.07465478  1.000000000 -0.003388206
## bisiesto     0.01410159 -0.05601013  0.02981413 -0.003388206  1.000000000

La matriz de correlación no muestra que exista una correlación especialmente alta entre ningún par de parámetros. Analicemos si se continúa cumpliendo la hipótesis de ruido blanco.

Box.test.2(residuals(fit.calendar),
           nlag = c(6, 12, 18, 24, 30, 36, 48))

##      Retard   p-value
## [1,]      6 0.5127498
## [2,]     12 0.4139199
## [3,]     18 0.3729273
## [4,]     24 0.3281915
## [5,]     30 0.5776594
## [6,]     36 0.3550016
## [7,]     48 0.6230107

Los p-values, aunque han disminuido, se han mantenido altos (mayores que 0,05) por lo que existe ruido blanco. Visualicemos mediante las gráficas de autocorrelación.

acf(fit.calendar$residuals, lag.max = 48, xlab = "Hipotecas",
     main = "Funcion de autocorrelacion simple")

Al observar la gráfica de autocorrelación simple se puede notar que sí existe ruido blanco tal como indicaban las correlaciones. En el análsis descriptivo inicial se pudo observar la existencia de datos atípicos, en la siguiente sección se determina de manera analítica si existe datos atípicos que puedan afectar la estabilidad del modelo en predicciones futuras.

Identificación de outliers

A continuación se determina analíticamente si existen casos atípicos en el modelo seleccionado.

outliers.train <- locate.outliers(fit.calendar$residuals,
                                      pars = coefs2poly(fit.calendar),
                                      types = c("AO", "LS", "TC"),cval=3)

outliers.train$abststat=abs(outliers.train$tstat)
outliers.train$abststat

## [1] 3.154992 3.168897 3.069912 3.933810 3.025998 3.008291

En este caso, como fue visto en clase, se buscan tres tipos de atipicidad. Los outliers de tipo AO son comparables a un impulso y afectan a la serie solamente en un momento puntual. Los outliers de tipo LS implican un cambio de nivel en la serie y los TC producen un cambio transitorio. Se ha identificado que el modelo es sensible a 6 casos atípicos. Asociemos los outliers identificados con el mes que le corresponden para facilitar la interpretación.

data.train$ind <- as.numeric(row.names(data.train))
outliers.train.date <- merge(outliers.train, data.train[,c('ind', 'date')],  by = "ind")
arrange(outliers.train.date,desc(outliers.train.date$abststat))

##   ind type   coefhat     tstat abststat       date
## 1 102   LS -31.21661 -3.933810 3.933810 2011-06-01
## 2  41   AO  27.98171  3.168897 3.168897 2006-05-01
## 3  26   AO -27.85893 -3.154992 3.154992 2005-02-01
## 4 146   AO  27.10766  3.069912 3.069912 2015-02-01
## 5 124   LS -24.01271 -3.025998 3.025998 2013-04-01
## 6 103   TC -26.17466 -3.008291 3.008291 2011-07-01

El outlier más significativo está asociado a un tipo LS (cambio de nivel) que se produjo en junio de 2011. Los siguientes más importantes son impulsos puntuales en mayo de 2006, febrero de 2005 y febrero de 2015. Seguidos de otro cambio de nivel en abril de 2013 y un cambio transitorio en julio de 2011.

Según la Asociación Hipotecaria Española (AHE), en el artículo “Las hipotecas registraron en el 2011 la mayor caída de la historia”, al cierre de 2011 se registra un desplome del saldo total de crédito hipotecario de las entidades financieras que fue considerado, en aquel momento, la mayor caída de la historia registrada. Este hecho pudiera asociarse al valor atípico observado en junio de 2011.

Asimismo, en 2006 La Caixa publica su Informe Mensual nº290, donde descarta que haya burbuja y vaticina una “desaceleración suave” del sector inmobiliario, lo cual podría asociarse a los outliers de 2006 y 2005.

Según el informe de la Sociedad de Tasación, en el año 2015 el precio de vivienda nueva permaneció estable en Comunidades como Aragón, Castilla-La Mancha, Castilla y León, Murcia y País Vasco, lo cual no se ajusta al valor atípico de febrero de 2015. Tampoco se ha encontrado una interpetación evidente para el resto de los casos por lo cual no serán incluidos en el modelo.

Modelo 6: ARIMA(0,1,1)(1,0,0)[12] + calendario + outliers

outliers.train <- outliers(c("LS","AO", "AO"), c(102, 41, 26))

outliers.variables.train <- outliers.effects(outliers.train, length(fit.calendar$residuals))
outliers.calendar.train <- as.matrix(cbind(calendar.train, outliers.variables.train))

fit.calendar.outliers <- Arima(data.train.ts,
                               order = c(0, 1, 1),
                               seasonal = list(order = c(1, 0, 0), period=12),
                               xreg = outliers.calendar.train,
                               method = 'ML')
fit.calendar.outliers

## Series: data.train.ts 
## Regression with ARIMA(0,1,1)(1,0,0)[12] errors 
## 
## Coefficients:
##           ma1    sar1  semanaSanta      dt  bisiesto     LS102     AO41
##       -0.6163  0.5369      -2.3897  1.3371   -0.4039  -30.8277  28.4323
## s.e.   0.0658  0.0599       0.5251  0.2476    4.3371    7.3879   8.4528
##           AO26
##       -30.5488
## s.e.    8.6516
## 
## sigma^2 estimated as 117.4:  log likelihood=-724.31
## AIC=1466.62   AICc=1467.62   BIC=1495.89

coeftest(fit.calendar.outliers)

## 
## z test of coefficients:
## 
##               Estimate Std. Error z value  Pr(>|z|)    
## ma1          -0.616346   0.065776 -9.3704 < 2.2e-16 ***
## sar1          0.536914   0.059855  8.9702 < 2.2e-16 ***
## semanaSanta  -2.389744   0.525075 -4.5512 5.333e-06 ***
## dt            1.337119   0.247631  5.3996 6.677e-08 ***
## bisiesto     -0.403878   4.337069 -0.0931 0.9258065    
## LS102       -30.827726   7.387902 -4.1727 3.010e-05 ***
## AO41         28.432311   8.452817  3.3636 0.0007692 ***
## AO26        -30.548773   8.651579 -3.5310 0.0004140 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Los coeficientes del modelo indican que las variables relativas a los outliers incluidas son significativas para el modelo, por lo tanto serán conservadas en el modelo final. Analicemos si su inclusión afectó la correlación entre los parámetros.

cor.arma(fit.calendar.outliers)

##                     ma1        sar1  semanaSanta          dt     bisiesto
## ma1          1.00000000  0.06978565 -0.054586078  0.02125063 -0.029603893
## sar1         0.06978565  1.00000000  0.014328896 -0.03354399 -0.117034394
## semanaSanta -0.05458608  0.01432890  1.000000000  0.07775538  0.030533877
## dt           0.02125063 -0.03354399  0.077755376  1.00000000 -0.011692521
## bisiesto    -0.02960389 -0.11703439  0.030533877 -0.01169252  1.000000000
## LS102        0.04615577  0.05237348  0.010377490 -0.01465218 -0.032724965
## AO41         0.08585830 -0.01154588  0.007822014 -0.06333532 -0.006172089
## AO26        -0.03224553 -0.10396188  0.031354211 -0.03474456  0.220635748
##                    LS102         AO41        AO26
## ma1          0.046155773  0.085858297 -0.03224553
## sar1         0.052373483 -0.011545884 -0.10396188
## semanaSanta  0.010377490  0.007822014  0.03135421
## dt          -0.014652182 -0.063335320 -0.03474456
## bisiesto    -0.032724965 -0.006172089  0.22063575
## LS102        1.000000000  0.004445407 -0.01104467
## AO41         0.004445407  1.000000000 -0.02845544
## AO26        -0.011044667 -0.028455439  1.00000000

Al analizar los coeficientes de correlación se evidencia que no existe un alto valor (mayor que 0.8 en valor absoluto) entre ninguno de los parámetros. Por lo tanto, no hay riesgo de colinealidad. Analicemos si se preserva la existencia de ruido blanco.

Box.test.2(residuals(fit.calendar.outliers),
           nlag = c(6, 12, 18, 24, 30, 36, 48))

##      Retard   p-value
## [1,]      6 0.6021314
## [2,]     12 0.6431322
## [3,]     18 0.7356284
## [4,]     24 0.9019731
## [5,]     30 0.9509610
## [6,]     36 0.8720311
## [7,]     48 0.9340538

acf(fit.calendar.outliers$residuals, lag.max = 48, xlab = "Hipotecas",
     main = "Funcion de autocorrelacion simple")

Tanto el test de Portemanteau como el gráfico de aucorrelación simple muestra que la serie es ruido blanco. Por lo tanto, este es nuestro modelo final. En la siguiente sección se evaluarán algunos de los modelos propuestos en función de las predicciones obtenidas en el conjunto de test.

Predicciones

Adicionamos las variables explicativas aplicadas anteriormente al conjunto de entrenamiento al conjunto de test.

exp.calendar.test <- calculoExplicativasCalendario(data.test$date,domingoYFestivosJuntos=0)

outliers.variables.test <- outliers.variables.train[181:192,]

outliers.calendar.test <- 
  as.matrix(cbind(exp.calendar.test[,c("semanaSanta", "dt", "bisiesto")], outliers.variables.test))

preds <- as.data.frame(predict(fit.calendar.outliers, n.ahead=12, newxreg=outliers.calendar.test))

Límites de confianza al 95%

U <- (preds$pred + 2*preds$se)**(1/lambda)
L <- (preds$pred - 2*preds$se)**(1/lambda)

data.pred <- data.frame(date = data.test$date, pred = (preds$pred)**(1/lambda),
                         LimSup = U, LimInf =L)
data.real <- merge(data[,c("date","total")], data.pred, by = "date", all.x = T)

Gráfico Real vs. Predicción

En el siguiente gráfico se dibujan los valores reales y las predicciones sobre el conjunto de test, que corresponde al último período de la serie que es el año 2019.

grafico1 <- ggplot(data = data.real) +
  geom_line(aes(x= date, y = total), color = 'steelblue',
            alpha = 0.8, size = 1) +
  geom_line(aes(x= date, y = pred), color = 'darkred',
            alpha = 0.9, linetype = 2, size = 1) + 
  xlab('FECHA') + ylab('Hipotecas')

ggplotly(grafico1)

Gráfico Real + Predicción + Límites

En el siguiente gráfico se dibujan los valores reales, las predicciones y los intervalos de confianza calculados anteriormente teniendo en cuenta la desviación estándar.

data.real$pred[year(data.real$date) != 2019] <- NA

grafico2 <- ggplot(data = data.real) +
  geom_line(aes(x= date, y = total), color = 'steelblue',
            alpha = 0.8, size = 0.8) +
  geom_line(aes(x= date, y = pred), color = 'darkred',
            size = 1)   +
  geom_line(aes(x= date, y = LimSup), color = 'orange',
            size = 1)  +
  geom_line(aes(x= date, y = LimInf), color = 'orange',
            size = 1) +
  xlab('FECHA') + ylab('Hipotecas')

grafico2

En la siguiente sección se comparan las predicciones que se obtendrían ajustando automáticamente un modelo a cada una de las provincias que integran la comunidad autónoma y sumando las predicciones que cada uno de ellos genera.

Ajuste automático

En el proceso automático se ajusta una serie para cada una de las provincias del País Vasco y posteriormente, se suman las predicciones individuales para obtener valores comparables con el resto de modelos desarrollados.

Provincia Alava

data.alava.train <- subset(original, year(original$date) != 2019)[c('date', 'Alava')]
data.alava.test <- subset(original, year(original$date) == 2019)[c('date', 'Alava')]

data.alava.train$exp <- data.alava.train$Alava**lambda
data.alava.test$exp <- data.alava.test$Alava**lambda

data.alava.train.ts <- as.ts(data.alava.train$exp)
data.alava.test.ts <- as.ts(data.alava.test$exp)

exp.calendar.alava.train <- calculoExplicativasCalendario(data.alava.train$date,domingoYFestivosJuntos=0)
exp.calendar.alava.test <- calculoExplicativasCalendario(data.alava.test$date, domingoYFestivosJuntos = 0)

calendar.alava.train <- 
  as.matrix(
    exp.calendar.alava.train[,c("semanaSanta", "dt", "bisiesto")])

calendar.alava.test <- 
  as.matrix(
    exp.calendar.alava.test[,c("semanaSanta", "dt", "bisiesto")])


alava.fit <- auto.arima(data.alava.train.ts,
                               max.d=1, max.D=1,
                               max.p=2, max.P=2,
                               max.q=2, max.Q=2, 
                               seasonal=TRUE,
                               ic="aic",
                               allowdrift=FALSE,
                               xreg=calendar.alava.train,
                               stepwise=TRUE)


alava.fit

## Series: data.alava.train.ts 
## Regression with ARIMA(0,1,1) errors 
## 
## Coefficients:
##           ma1  semanaSanta      dt  bisiesto
##       -0.7496      -1.1844  0.1075    3.3785
## s.e.   0.0558       0.5220  0.2531    4.9882
## 
## sigma^2 estimated as 115.7:  log likelihood=-723.15
## AIC=1456.31   AICc=1456.63   BIC=1472.57

cor.arma(alava.fit)

##                       ma1 semanaSanta            dt     bisiesto
## ma1          1.0000000000 -0.13731307 -0.0002557799  0.004602419
## semanaSanta -0.1373130694  1.00000000  0.0309148878  0.062263884
## dt          -0.0002557799  0.03091489  1.0000000000 -0.004746390
## bisiesto     0.0046024191  0.06226388 -0.0047463898  1.000000000

Box.test.2(residuals(alava.fit),
           nlag = c(6, 12, 18, 24, 30, 36, 48))

##      Retard    p-value
## [1,]      6 0.00518674
## [2,]     12 0.00005911
## [3,]     18 0.00000236
## [4,]     24 0.00000005
## [5,]     30 0.00000053
## [6,]     36 0.00000007
## [7,]     48 0.00000000

acf(alava.fit$residuals, lag.max = 48, xlab = "Hipotecas",
     main = "Funcion de autocorrelacion simple")

alava.preds <- as.data.frame(predict(alava.fit, n.ahead=12, newxreg=calendar.alava.test))

Provincia Gipuzkoa

data.gipuzkoa.train <- subset(original, year(original$date) != 2019)[c('date', 'Gipuzkoa')]
data.gipuzkoa.test <- subset(original, year(original$date) == 2019)

data.gipuzkoa.train$exp <- data.gipuzkoa.train$Gipuzkoa**lambda 
data.gipuzkoa.test$exp <- data.gipuzkoa.test$Gipuzkoa**lambda

data.gipuzkoa.train.ts <- as.ts(data.gipuzkoa.train$exp, frequency = 12)
data.gipuzkoa.test.ts <- as.ts(data.gipuzkoa.test$exp, frequency = 12)

exp.calendar.gipuzkoa.train <- calculoExplicativasCalendario(data.gipuzkoa.train$date,domingoYFestivosJuntos=0)
exp.calendar.gipuzkoa.test <- calculoExplicativasCalendario(data.gipuzkoa.test$date, domingoYFestivosJuntos = 0)

calendar.gipuzkoa.train <- 
  as.matrix(
    exp.calendar.gipuzkoa.train[,c("semanaSanta", "dt", "bisiesto")])

calendar.gipuzkoa.test <- 
  as.matrix(
    exp.calendar.gipuzkoa.test[,c("semanaSanta", "dt", "bisiesto")])


gipuzkoa.fit <- auto.arima(data.gipuzkoa.train.ts,
                               max.d=1, max.D=1,
                               max.p=2, max.P=2,
                               max.q=2, max.Q=2, 
                               seasonal=TRUE,
                               ic="aic",
                               allowdrift=FALSE,
                               xreg=calendar.gipuzkoa.train,
                               stepwise=TRUE)

gipuzkoa.fit

## Series: data.gipuzkoa.train.ts 
## Regression with ARIMA(0,1,2) errors 
## 
## Coefficients:
##           ma1      ma2  semanaSanta      dt  bisiesto
##       -0.6153  -0.1240      -0.7658  0.6168    6.2268
## s.e.   0.0728   0.0696       0.3783  0.1736    3.6619
## 
## sigma^2 estimated as 63.56:  log likelihood=-665.35
## AIC=1342.71   AICc=1343.17   BIC=1362.22

cor.arma(gipuzkoa.fit)

##                     ma1         ma2 semanaSanta          dt    bisiesto
## ma1          1.00000000 -0.73477444 -0.12083657  0.05472761 -0.15642898
## ma2         -0.73477444  1.00000000  0.08932132 -0.06401154  0.14489691
## semanaSanta -0.12083657  0.08932132  1.00000000  0.02190132  0.09776195
## dt           0.05472761 -0.06401154  0.02190132  1.00000000 -0.01605143
## bisiesto    -0.15642898  0.14489691  0.09776195 -0.01605143  1.00000000

Box.test.2(residuals(gipuzkoa.fit),
           nlag = c(6, 12, 18, 24, 30, 36, 48))

##      Retard    p-value
## [1,]      6 0.90063518
## [2,]     12 0.00037983
## [3,]     18 0.00030865
## [4,]     24 0.00001086
## [5,]     30 0.00001261
## [6,]     36 0.00000786
## [7,]     48 0.00000636

acf(gipuzkoa.fit$residuals, lag.max = 48, xlab = "Hipotecas",
     main = "Funcion de autocorrelacion simple")

gipuzkoa.preds <- as.data.frame(predict(gipuzkoa.fit, n.ahead=12, newxreg=calendar.gipuzkoa.test))

Provincia Vizcaya

data.vizcaya.train <- subset(original, year(original$date) != 2019)[c('date', 'Vizcaya')]
data.vizcaya.test <- subset(original, year(original$date) == 2019)[c('date', 'Vizcaya')]

data.vizcaya.train$exp <- data.vizcaya.train$Vizcaya**lambda 
data.vizcaya.test$exp <- data.vizcaya.test$Vizcaya**lambda

data.vizcaya.train.ts <- as.ts(data.vizcaya.train$exp, frequency = 12)
data.vizcaya.test.ts <- as.ts(data.vizcaya.test$exp, frequency = 12)

exp.calendar.vizcaya.train <- calculoExplicativasCalendario(data.vizcaya.train$date,domingoYFestivosJuntos=0)
exp.calendar.vizcaya.test <- calculoExplicativasCalendario(data.vizcaya.test$date, domingoYFestivosJuntos = 0)

calendar.vizcaya.train <- 
  as.matrix(
    exp.calendar.vizcaya.train[,c("semanaSanta", "dt", "bisiesto")])

calendar.vizcaya.test <- 
  as.matrix(
    exp.calendar.vizcaya.test[,c("semanaSanta", "dt", "bisiesto")])


vizcaya.fit <- auto.arima(data.vizcaya.train.ts,
                               max.d=1, max.D=1,
                               max.p=2, max.P=2,
                               max.q=2, max.Q=2, 
                               seasonal=TRUE,
                               ic="aic",
                               allowdrift=FALSE,
                               xreg=calendar.vizcaya.train,
                               stepwise=TRUE)

vizcaya.fit

## Series: data.vizcaya.train.ts 
## Regression with ARIMA(0,1,1) errors 
## 
## Coefficients:
##           ma1  semanaSanta      dt  bisiesto
##       -0.6111      -1.4653  1.3248   -0.6452
## s.e.   0.0648       0.4362  0.2039    4.1663
## 
## sigma^2 estimated as 87.39:  log likelihood=-696.15
## AIC=1402.31   AICc=1402.63   BIC=1418.57

cor.arma(vizcaya.fit)

##                      ma1 semanaSanta           dt     bisiesto
## ma1          1.000000000 -0.16316380 -0.004351635 -0.023090197
## semanaSanta -0.163163801  1.00000000  0.032229442  0.084767888
## dt          -0.004351635  0.03222944  1.000000000 -0.003801747
## bisiesto    -0.023090197  0.08476789 -0.003801747  1.000000000

Box.test.2(residuals(vizcaya.fit),
           nlag = c(6, 12, 18, 24, 30, 36, 48))

##      Retard    p-value
## [1,]      6 0.20941483
## [2,]     12 0.00008122
## [3,]     18 0.00011075
## [4,]     24 0.00001618
## [5,]     30 0.00001367
## [6,]     36 0.00000721
## [7,]     48 0.00002919

acf(vizcaya.fit$residuals, lag.max = 48, xlab = "Hipotecas",
     main = "Funcion de autocorrelacion simple")

viscaya.preds <- as.data.frame(predict(vizcaya.fit, n.ahead=12, newxreg=calendar.vizcaya.test))

Los modelos ajustados automáticamente para las provincias Alava, Gipuzkoa y Vizcaya fueron ARIMA(0,1,1), ARIMA(0,1,2) y ARIMA(0,1,1) respectivamente, en ninguno de los casos se incluye la parte estacional. El modelo de la provincia de Alava fue ARIMA(0,1,1), sin embargo, tanto el test como la gráfica de autocorrelación simple muestra que no cumple la hipótesis de ruido blanco. Así mismo ocurre con las provincias de Gipuzkoa y Vizcaya. Al observar los gráficos de autocorrelación simple se puede observar en todos los casos que los palitos que más sobresalen son múltiplos de 12, por lo tanto, incluir una parte estacional con período 12, como se hizo en el modelo propuesto en este trabajao, podría llevar a la existencia de ruido blanco en estos casos. Para la construcción de un modelo comparable al modelo final construido se suman las predicciones de cada provincia de la comunidad.

data.pred.alava <- data.frame(date = data.test$date, pred = (alava.preds$pred)**(1/lambda))
data.pred.gipuzkoa <- data.frame(date = data.test$date, pred = (gipuzkoa.preds$pred)**(1/lambda))
data.pred.viscaya <- data.frame(date = data.test$date, pred = (viscaya.preds$pred)**(1/lambda))

data.pred.sum <- data.frame(date = data.test$date, pred_sum = (data.pred.alava$pred + data.pred.gipuzkoa$pred + data.pred.viscaya$pred)) 
all.preds <- merge(data.real, data.pred.sum, by = "date", all.x = T)

En el siguente gráfico se presentan los valores reales (color azul), predicciones (color rojo) y las predicciones resultantes de la suma (línea verde).

grafico3 <- ggplot(data = all.preds) +
  geom_line(aes(x= date, y = total), color = 'steelblue',
            alpha = 0.8, size = 1) +
  geom_line(aes(x= date, y = pred), color = 'darkred',
            alpha = 0.9, linetype = 2, size = 1) + 
  geom_line(aes(x= date, y=pred_sum), color = '#52854C',
            alpha = 0.9, linetype = 2, size = 1) +
  xlab('FECHA') + ylab('Hipotecas')

ggplotly(grafico3)

Visualmente se puede notar que la curva correspondiente al modelo ajustado ARIMA se aproxima más a los valores reales que el obtenido mediante auto.arima en los primeros meses de 2019, sin embargo, en los últimos meses nuestro primer modelo se aleja más que el modelo automático. Comprobemos mediante el cálculo analítico de los errores.

Cálculo de errores

En esta sección se calcula el error de tres de los modelos presentados:

• ARIMA(0, 1, 1)x(1,0,0)[12]
• ARIMA(0, 1, 1)x(1,0,0)[12] + calendario + outliers
• AUTOARIMA

Error in-sampling

preds.arima <- as.data.frame(fit.4$fitted)**(1/lambda)
names(preds.arima) <- c("ARIMA(0, 1, 1)x(1,0,0)[12]")

preds.calendar.outliers <- as.data.frame(fit.calendar.outliers$fitted)**(1/lambda)
names(preds.calendar.outliers) <- c("ARIMA(0, 1, 1)x(1,0,0)[12] + calendario + outliers")

preds.autoarima <- as.data.frame(alava.fit$fitted)**(1/lambda) + as.data.frame(gipuzkoa.fit$fitted)**(1/lambda) + as.data.frame(vizcaya.fit$fitted)**(1/lambda)
names(preds.autoarima) <- c("AUTOARIMA")

real_predictions <- cbind(data.train,preds.arima, preds.calendar.outliers, preds.autoarima)

real_predictions$monthly.MAPE.arima <- abs(100*(real_predictions$total - 
                                                   real_predictions$`ARIMA(0, 1, 1)x(1,0,0)[12]`)/real_predictions$total)


real_predictions$monthly.MAPE.calendar.outliers <- abs(100*(real_predictions$total - 
                                                   real_predictions$`ARIMA(0, 1, 1)x(1,0,0)[12] + calendario + outliers`)/real_predictions$total)


real_predictions$monthly.MAPE.autoarima <- abs(100*(real_predictions$total - 
                                                   real_predictions$AUTOARIMA)/real_predictions$total)

Error global in-sampling sobre histórico:

global.MAPE.arima <- mean(real_predictions$monthly.MAPE.arima)
print(paste0("ARIMA(0, 1, 1)x(1,0,0)[12]: ", global.MAPE.arima))

## [1] "ARIMA(0, 1, 1)x(1,0,0)[12]: 11.9665242739954"

global.MAPE.calendar <- mean(real_predictions$monthly.MAPE.calendar.outliers)
print(paste0("ARIMA(0, 1, 1)x(1,0,0)[12] + calendario + outliers: ", global.MAPE.calendar))

## [1] "ARIMA(0, 1, 1)x(1,0,0)[12] + calendario + outliers: 10.4870722168623"

globa.MAPE.autoarima <- mean(real_predictions$monthly.MAPE.autoarima)
print(paste0("AUTOARIMA: ", globa.MAPE.autoarima))

## [1] "AUTOARIMA: 12.820226544532"

Error out-of-sampling

Calculemos los errores sobre el conjunto de test

preds <- all.preds[193:204,]
preds_fit.4 <- as.data.frame(predict(fit.4, n.ahead=12)$pred)**(1/lambda)
colnames(preds_fit.4) <- c("pred.arima")

preds <- cbind(preds, preds_fit.4)
colnames(preds) <- c("date", "real", "calendar.outliers.pred", "LimSup", "LimInf", "autoarima.pred", "arima.pred")

preds$monthly.MAPE.arima <- abs(100*(preds$real - preds$arima.pred)/preds$real)
preds$monthly.MAPE.calendar.outliers <- abs(100*(preds$real - preds$calendar.outliers.pred)/preds$real)
preds$monthly.MAPE.autoarima <-  abs(100*(preds$real - preds$autoarima.pred)/preds$real)

Error out-of-samplig sobre test.

test.MAPE.arima <- mean(preds$monthly.MAPE.arima)
print(paste0("ARIMA(0, 1, 1)x(1,0,0)[12]: ", test.MAPE.arima))

## [1] "ARIMA(0, 1, 1)x(1,0,0)[12]: 9.25247746645789"

test.MAPE.calendar.outliers <- mean(preds$monthly.MAPE.calendar.outliers)
print(paste0("ARIMA(0, 1, 1)x(1,0,0)[12] + calendario + outliers: ", test.MAPE.calendar.outliers))

## [1] "ARIMA(0, 1, 1)x(1,0,0)[12] + calendario + outliers: 8.62845721290206"

test.MAPE.autoarima <- mean(preds$monthly.MAPE.autoarima)
print(paste0("AUTOARIMA: ", test.MAPE.autoarima))

## [1] "AUTOARIMA: 9.86055619957073"

Al analizar los errores se puede observar que el modelo que menor error alcanza es el ARIMA(0, 1, 1)x(1,0,0)[12] incluyendo las varibles explicativas de calendario y valores atípicos, tanto en el conjunto de test como en el conjunto de training. Le sigue el modelo ARIMA(0,1,1)x(1,0,0)[12] sin incluir el calendario y valores atípicos y por último, el modelo de suma de las provincias. Sin embargo, la diferencia entre errores es muy pequeña.

Conclusiones

En este trabajo se han desarrollado un conjunto de modelos ARMA para el análisis de datos temporales. La selección de modelos se ha hecho siguiendo la metodología Box-Jenkins estudiada en clase. Tras el , se puede concluir que el modelo que mejor resultados alcanzó fue el ARIMA(0, 1, 1)x(1,0,0)[12] que incluye las varibles explicativas de calendario y valores atípicos, aunque la diferencia respecto al resto de los modelos no es muy diferenciadora. Esto pudiera sugerir que la inclusión de variables explicativas, que aumentan la complejidad del modelo y el tiempo invertido, no hacen mucha diferencia. Sin embargo, debemos tener en cuenta que se está midiendo la efectividad de los modelos en un conjunto con solo 12 muestras correspondientes a un solo año. Por lo cual la inclusión de variables explicativas es recomendable, porque permiten simular mejor el comportamiento de datos futuros al anadir información externa. Por otra parte, el enfoque basado en la suma de predicciones mostró un rendimiento más bajo aún incluyendo variables relativas al calendario lo cual sugiere, en este caso, que un análisis personalizado y adaptado al escenario particular arroja mejores resultados que los procesos auomáticos.